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HOW TO MODEL THE COVARIANCE STRUCTURE IN A 
SPATIAL FRAMEWORK: VARIOGRAM OR CORRELATION 

FUNCTION? 

GIOVANNI PISTONE AND GRAZIA VIGARIO 


Abstract. The basic Kriging’s model assumes a Gaussian distribution 
with stationary mean and stationary variance. In such a setting, the joint 
distribution of the spatial process is characterized by the common variance 
and the correlation matrix or, equivalently, by the common variance and 
the variogram matrix. We discuss in in detail the option to actually use 
the variogram as a parameterization. 

Kriging, Universal Kriging, Non-parametric Kriging, Elliptope, Bayes 

1. Introduction 

The present note is a development of the Author’s paper [7]. An application 
was implemented in Vicario et al. paper [12]. In the interest of clarity we allow 
here for a little overlap with the papers referred to above. 

We discuss the notion of variogram as it is used in Geostatistics and we offer 
some preliminary thought about the possibility of a non-parametric approach 
to Universal Kriging aiming to the use of a Bayes methodology. Variograms 
are due to Matheron [6], and there are many modern expositions i.e., Cressie 
[3, Ch. 2], Chiles and Delhner [2, Ch. 2], Gneiting et al. [5], Gaetan and 
Guyon [4, Ch.lj. 

In Sec. 2 we give a brief overview of the so called Universal Kriging model 
and its parameterization with the Matheron’s variogram function. In Sec. 3 
we dehne the general variogram matrix and give a necessary and sufficient 
condition for a positive variance and a matrix T to be a variogram matrix 
of a covariance a^R, where i? is a correlation matrix. In Sec. 4 we provide 
some useful computations concerning the inverse variogram matrix. Sec 5 is 
devoted to an interpretation of the variogram matrix as related to a projection 
of the Gaussian held. In Sec. 6 we discuss the shape of the set of parameters 
of the general Kriging model. A section of conclusions closes the paper. 

2. Universal Krige setup 

We consider a Gaussian n-vector K, n > 2, whose mean has the form 
1-1 = fil and whose covariance matrix S = has constant diagonal 

* = 1, • • •, R- The assumption on the mean and the diagonal terms 
is the weakest stationarity assumption, i.e. a 1-st order stationarity. We can 
write Y ~ N„(/il, where p is a general mean value, the common 

variance, and R = [pij]ij=i a generic correlation matrix. 
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The variogram of Y is the n x n matrix T = ['Jij]^j=i whose element jij 
is half the variance of the difference Yi — Yj. As the mean value is constant, 
the variance of the difference is equal to the second moment of the difference. 
It is expressed in terms of the common variance and the correlations pij, 

j = 1 ,... ,n, as 

2'yij = Var (Yi - Yj) = 0 -^( 6 ^ - ej)'R{ei - Sj) = 

^ {.Pa T Pjj ‘^Pij) (1 Pij) ; 

and, in matrix form, as 

T = a‘^{ll' - R). 

The simple Gaussian model described above is commonly used in Geostatis¬ 
tics, when each random component Y^ of the random vector Y is associated 
to a location Xj, i = 1,..., n, in a given region X, G X, i = 1,..., n. 

We briefly describe the most common set-up in Geostatistics. The elements 
of the variogram matrix 7 ^ are assumed to be a given function 7 of the distance 
between two locations, 7 *^ = 'y{d{xi,Xj)). In such a case the statistical model 
is characterized by the choice of a distance d{x, y), x,y E X, and by the choice 
of a function 7 , called variogram function, dehned on a real domain containing 
{d{x,y)\x,y G X} . The existence of a positive and a correlation matrix 
R = [pij] such that < 7^(1 — pij) = 'y{d{xi,Xj)) imposes a nontrivial condition 
on the function 7 , see Sasvari [11] and Gneiting et al. [5] and []. Such a 
model, where it is assumed that the vector of means is constant, pL = p,! is 
called universal Krige model. We do not consider the more general case of a 
non-constant mean. 

Krige has further qualihed this model by adding assumptions on the vari¬ 
ance function and suggesting a statistical method to estimate the value at 
an untried point xq given a set of observation Yi,... ,Yn at points xi,..., Xn- 
Precisely: 

(1) Krige’s modeling idea is to assume the variogram function 7 to be 
an increasing function on [ 0 , cx)[, so that the variogram’s values are 
increasing with the distance. Moreover, the correlation between loca¬ 
tions is assumed to be positive. The rational is to model a variability 
which increases with the distance and is bounded by a general variance: 

0 < ^ Var {Yi - Yj) = jij = 'y{d{xi, Xj)) =a‘^{l- pij) < . 

The increasing function 7 : M_|_ —)■ M_|_ is assumed to be continuous 
everywhere but at 0 . As it is bounded at -|-cx), the general shape is as 

in Fig. 1. 

The parameters in the Krige’s universal model are unrestricted val¬ 
ues of /i G M, > 0, and restricted values of R that are usually 
estimated over a suitable parametric model. 

(2) Krige’s idea is to predict the value Yq = Y^q at an untried location 
Xq with the conditional expectation based on a plug-in estimate of 
the parameters. If I = {1,... ,n} and the locations in the model are 
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Figure 1. A general variogram function is 0 at 0, can have 
jump at 0 which is called Nugget, has a hnite limit at +oo 
named Sill. The Range is a length such that the value is equal 
to the limit value for any practical purpose. 

Xo,Xi,..., Xn, the regression value is 

= with S= ^ 2 ^’° 

The set of data that give the same prediction is an affine plane in M”. 
The variance of the prediction is (Jq — 

In this paper we do not follow this approach, but we adopt a general non- 
parametric attitude, where /i is real number, is a positive real number, R 
is a positive dehnite matrix with unit diagonal, possibly with positive entries. 
The variogram matrix T is not restricted and we do not enforce the existence 
of any special form. 


3. The variogram matrix 

Our plan now is to express the Krige’s computations in term of Matheron’s 
variogram matrix T. One good reason to use T as a basic parameter is because 
its empirical estimator is unbiased. 

Let us discuss in some detail the basic transformation of matrix parameters 
(1) T = - R) = aHl'-Y . 

We note that T = 0 if, and only if, R = 11', such extreme case being always 
excluded in the following. In fact, in most cases we will assume deti? ^ 0. 

The entries of T are nonnegative and bounded by 2a^ because the correla¬ 
tions are bounded between —1 and 1. If all the correlations are nonnegative, 
the the entries of T are bounded by a^. 
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The difference between the covariance matrix and the variogram matrix is 
a matrix of rank 1, T + S = Let us remark moreover that 

( 2 ) iii' = _L(r + E) 

n na^ 

is the orthogonal projector on the space of constant vectors Span (1). 

We denote by S= the cone of nonnegative dehnite matrices with constant 
diagonal, by S=i the convex set of correlation matrices, and by V the cone of 
variogram matrices. We have the following characterization of V. 

Proposition 1. A nonzero matrix T is the variogram matrix of some covari- 
anee matrix of the form S = a^R, with > 0 and R a correlation matrix, if, 
and only if, the three following conditions hold 

( 1 ) r is symmetric, and has zero diagonal; 

(2) r is conditionally negative dehnite, i.e. wTw < 0 if w'l = 0 ; 

(3) sup {a;Ta;|a;'l = 1} < 

Proof. Assume T = — i?) 7 ^ 0 with R correlation matrix and > 0. 

Condition 1 follows from the dehnition. If we write a generic vector as a; = 
w + al with w'l = 0 , we have 

= x'Tx + x'Jlx . 

In particular. Condition 2 follows because a = 0 implies a^Ta; = —x'T^x. 
Finally, if x'l = 1, that is na = 1, we have x'Tx — = x'Hx > 0 and 

Condition 3 follows. 

Viceversa, consider the matrix 11' — a~^r. It is symmetric, with unit 
diagonal. We need only to show it is positive dehnite: 

a;'(ll' — a~^r)x = {x'lY — a~^x'Tx = 

because of Condition 3. □ 


The lower bound imposed on means that the parameterization with 
carrying one degree of freedom, and with T, carrying n{n — l )/2 degrees of 
freedom, has a drawback in that the two parameters are not independently 
dehned on a product set. Note the relation I'Tl = a‘^{n^ — VRT) that we are 
going to discuss in Sec. 4 below. 

In conclusion, there is a one-to-one transformation of parameters 

(a^T) fA AA S 


with G M>, T G V, i? G S=i, S G S=, namely: 

(1) The mapping from S G S= to the couple (a^, T) G M> x V factors as 


n 


§= 3 S ^ (- Tr (S), - Tr (S) S) = (a^ R) g] 0, cx)[x7^ 


n 


-2 
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and 

]0, oo[xTZ3 (a^, R) HA (a^, <7^(11' — R)) = 

(a^r) e {(a^r)|r e V,sup{a;'ra;|a;'l = 1} < (T^ . 

(2) Inverse is 

e V,sup{a;'ra;|a;'l = 1} < a^} 3 (a^ T) ^ 

a^ll' - r = S G §= 

4. Inverse variogram matrix T~^ 

The equations two above, both based on the dehnition of the variogram 
matrix as T = a‘^ {11' — R), provide a simple connection between the param¬ 
eterization based on the covariance matrix S and the parameterization based 
on the couple and T. However, we want to spell out the computation of 
an other key statistical parameter, namely the concentration matrix We 
begin by recalling a well known equation in matrix algebra [ 8 ]. We review the 
result in detail as we need an exact statement of the conditions under which 
is true in our case. 

Proposition 2 (Sherman-Morrison formula). Assume the matrix A is invert¬ 
ible. The matrix 11' — H is invertible if, and only if, l'H“^l 7 ^ 1. /n such a 
case, 

det(lT -A) = (-l)'^(l - I'A-H) det A , 

(IT - A)-^ = -A-^ - (1 - TA-^l)-^A-^lTA-\ 

Proof. The multi-linear expansion of det (11' — A) is written in terms of the 
adjoints {—A)'^ of each element {—A)ij by 


det(ll' — A) = det (—H) -I- 

j=i i=i 

n 

= (-1)" det A - (-1)"-^ ^ A'^ 

Li=i 

= (-1)" det A - (-l)''-il(adj A)l' . 

As det A 7 ^ 0, we can factor-out (—l)"'(det A) to get 

det(lT - A) = (-l)'^(det A)(l - TA'^) 
and the statement about the determinant follows. The inversion formula is 
directly checked. □ 

We are concerned with the invertibility of T = {11' — R), hence we need 

to discuss the condition l'i?“^l 7 ^ 1 . 

Proposition 3. Let R be a correlation matrix and assume detR 7 ^ 0. Let 
\j > 0, j = 1,... ,n, be the spectrum of R and let Uj be a set of unit eigen¬ 
vectors. It can be proved: 
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(1) Tri? = Yl"j=i ^ det-R = YYj^ -_i Aj < 1, with equality if, and 
only if, R = I. 

(2) Tri?“^ = ^ equality if, and only if, R = I. 

(3) ^ 1. 

Proof. (1) n = TtR = ^j- From det R = ]Yj=i the arithmetic 

mean is larger than the the geometric mean, 

\ / n \ n 

l = =(detfi)i. 

with equality if, and only if the A/s are all equal, hence equal to 1, 
which happens if i? = /. 

(2) The geometric mean is larger or equal than the harmonic mean, hence 

( n \ n / n 

n j ^ ^ 

with equality if, and only if, Aj = 1, j = 1 ,..., n. It follows ^ YTj=l ^ ^ 

1 . 

(3) We derive a contradiction from 1 = As R~^ = 

and E;=i(lS'f = l|l||' = r^^ 

n n 

1 = I'fl-a = ^at‘(iA)" = . 

i=i j=i 

where 6j = (l'w)^/n^ > 0 and 'YTj=i^j ~ 1- From the convexity of 
A H-)■ A“^ we obtain 

n ( 

j=l \j=l 

hence the contradiction 

1 . ^ 1 1 
^ < ;^max{A,|j = l,...,n} < - . 




□ 


From Proposition 3 we have immediately the following result of interest. 

Proposition 4. Assume the eorrelation matrix S = a^R G §= is invertible. 
It follows that r = (T^(ll' — R) is invertible, with 

(3) r-^ = - (a-2 - l'S-^l)-^S-ilTS-i 

and 

(4) = -F"^ - (ct-2 - TF-4)-^F-4TF-^ 
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Proof. From the assumption it follows det i? 7 ^ 0 so that 

^_2 _ ^ ^ Q ^ 

hence the conclusion. □ 

We can now analyze the likelihood of the Gaussian model N{fil,a^R) in 
terms of the variogram. 

First, we compute the determinant of the correlation matrix 
det = det (a^ll' — F) 

= det (ll'-a-^F) 

= [det (—(T“^F) + I'adj (—(T“^F) l] 

= det (—F) — adj (—F) 1 . 

Second, we compnte the qnadratic form of the concentration matrix 
= y' (-F"! - (a-^ - VT-Hr^T-Hl'T-^) y 

= -y'T-^y-{a-^-rT-H)-\y^T-Hr. 

Third, we compute the log-likelihood with p = 0: 

\ogp{y\a^,T) = 

77 1 1 

- 2 log ( 2 ^) “ 2 ~ ~ = 

77 1 

- 2 log ( 2 ^) “ 2 ~ 

+ ly'T-^y + - i'r-ii)-i(7/'F-ii)2. 

Here are the essentials of the computations leading to a maximnm likelihood 
estimation of F are the following. In the direction of a generic symmetric 
matrix with zero diagonal H, 

dniT ^ log (det (ll' - n-^F))) = Tr ((a^ll' - T)-^H) 

and 

dniT^yfir -a-^T)-^y) 

= Tr ((a^ll' - r)~^yy\a‘^ir - Vy^H) . 

so that the normal equations for F reduce to the condition 

— (a^ll' — F)“^ -I- (a^ll' — T)~^yy'{yil' — F)“^ is diagonal . 

The approach with parameters a^, F is feasible in principle, but it does not 
appear promising in term of ease of computation. We will see in Sec. 5 below 
a different, possibly better, approach. 
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5. Projecting on Span(l)‘^ 

We now change our point of view to consider the same problem from a 
different angle suggested by the observation the variogram does not change if 
we change the general mean /i. In fact, we can associate the variogram with 
the state space description of the Gaussian vector. 

The following proposition is a new characterization of the variogram matrix 
in our setting. 

Proposition 5. (1) The matrix T is a variogram matrix of a covariance 

matrix S G V+ if, and only if, the matrix 

(5) E„ = -(/-iii')'r(/-iir) 

is symmetric, positive definite and with constant diagonal. 

(2) If Yo N„(0,So), then its variogram is T and Yq is supported by 
Span (I)"*". 

Proof. (1) If r = {11' — R) is the variogram matrix of S = a^R, then 

from Eq. (5) we have 



which is indeed positive dehnite. Let us show that the diagonal ele¬ 
ments of So are constant. 

(So).. = a^e' - ^IT^ r(^I- ^1T^ e. 

/ 2 1 

= I e'iRei -e.i?l -|— 

\ n 

Viceversa, assume So is a covariance matrix. As Si — ej E Span (1)'*', 
the variogram of Sq has elements 

(Oj Cj'j So( 6 j Gj') 

(e. - e,y - in') (-r) - in') (e, - ej) = 

(o. Gj') r(ej 6 j) y.. 'yjj 2^.^. 

(2) As l'(e. — Gj) = 0 , then V (J — ^H^)^ (~r) (J — ^11') 1 = 0 , hence 
the distribution of Yq is supported by the space Span (I)"*". 

□ 
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It is possible to split every Gaussian Y with covariance matrix S G §= 
according the splitting M"" = Span(l) © Span (!)■*■. The corresponding pro¬ 
jections split the Gaussian process Y into two components, one with the co- 
variance as in Proposition 5, the other proportional to the empirical mean. 
Note that the two components have a singular covariance matrix. 


Proposition 6. Let Y ~ S), S = a^R G S+ with variogram V = 

— R). Let Y = (^I — ^11') F ~ (0, Sq) be the projection ofY onto 

Span (!)■'■ so that we can write Y = Y + Y, where each component ofY is the 
empirical mean ^I'Y. 

(1) The distribution ofY depends on the variogram only, 



and the variogram matrix ofY isT. 

(2) The distribution of ^'^'Y, conditionally to Y is Gaussian with mean p. 


Proof. 




□ 


This suggests the following empirical estimation algorithm. 

(1) Project the independent sample data yi,...,yp^ onto Span(l)''' by 

subtracting the empirical mean to get yi = yi — 

y+,... ,yN = y^ — y+. Use the empirical estimator of the variogram 
matrix on the projected data. 

(2) Estimate y with the empirical mean. 

In the same spirit, we could suggest the simulation of a random variable 
with variogram matrix P by the generation of N„(0, Sq) data in Span (I)"*". 

Both suggestions will be further discussed in future work. 


6. Elliptope 


We now turn to the geometrical description of the set of variograms. From 
the basic equation P = —R) it follows that the set of variogram matrices 

is an affine image of the set of correlation matrices in the space of symmetric 
matrices. The set of correlation matrices is a convex bounded set whose 
geometrical shape has been studied in a number of papers, i.e. [10], [10], [9]. 
Such a shape is of central interest in a non parametric Bayesian approach 
to the statistics of the universal Kriging model. It appears also in convex 
optimization, where it has been called elliptope. 

Let us discuss the case n = 3. All principal minors of R are nonnegative. 


/ 

'l 

X 

y 

det (R) = det I 

X 

1 

z 

V 

y 

2 : 

1 


1 — x'^ — y‘^ — z'^ + 2xyz > 0 
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1.0 


1,C 


0 



1.0 


Figure 2. The 
3-elliptope 


Figure 
3. Algebraic 
variety 


and 1 — — ip',1 — > 0. The three last inequalities dehne the cube 

Q = [— 1 , + 1 ]^ while the equation 1 — + 2xyz = 0 is a cubic 

algebraic variety whose intersection with the cube Q is the border of the 
elliptope. All horizontal section z = c, —l<c<l, of the elliptope are the 
interior of the ellipses 


1 — x"^ — y'^ + 2cxy > 


Same for other sections. See Fig.s 2 and 3. 

Various proposals of apriory distribution on the elliptope exist, see for ex¬ 
ample [ 1 ]. 

Going on with the discussion of our example, the volume is easily computed 
and the uniform apriori is defined. Simulation is feasible for example by 
rejection method. An other option is to write R = A' A where the columns of A 
are unit vectors. This gives an other possible apriori starting from independent 
unit vectors. Simulation is feasible, for example starting with independent 
standard Gaussians. 

An interesting option is the Gholesky representation. A symmetric matrix 
A is positive definite if there exists an upper triangular matrix 



T — 12 — 0 ^22 ^23 5 tii > 0 

tg 0 0 ^33 


such that 



Moreover, ^ 11 ^ 22^3 7 ^ 0 T is unique and invertible if, and only if, A is 
invertible. It is an identifiable parameterization for non singular matrices. 
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In the case of a the correlation matrix R = T’T with 



ti 


a/1 - Us 

T = 

t '2 

= 

0 ~ ^23 Us 


t's 


0 0 1 _ 


ti e 0* ^ X S, 


n-l+l 


It follows 


R = 


1 


vT 


V^l - 


'-12 


ffstis 1 

\/l ~ ^12 ~ ^13^13 ^12^13 + a/I — ^23^23 


\/l - 4 -4^1 3 

^12^13 + \/l ~ ^23^ 
1 


23 


and 


det (-R) — (1 ti 2 ^ 13 ) (1 ^ 23 ) 


7. Conclusion 

In the last decades Kriging models have been recommended not only for the 
original application, but spatial noisy data in general. Thanks to the availabil¬ 
ity of comprehensive computing facilities and the recent progresses in software 
development, the numerical simulation of technologically complex systems has 
become an attractive alternative option to the physical experimentation. The 
most popular meta-model used when dealing with Computer Experiments 
(CE) is the Kriging model. The accuracy of this model strongly depends on 
the detection of the correlation structure of the responses. In the Bayesian 
approach, where the posterior distribution of a prediction Krige’s To given the 
training set (Ti,..., K„) requires less uncertainty as possible on the correla¬ 
tion function, the use of the variogram as a parameter should be preferred 
because it does not demand a parametric approach as the correlation estima¬ 
tion does. The authors proved in a previous paper [7] the equivalence between 
variogram and spatial correlation function for stationary and intrinsically sta¬ 
tionary processes. Here the study has been devoted to the characterization 
of matrices which are admissible variograms in the case of 1-stationarity. We 
expect these Endings will allow for the imputation of an apriori distribution 
on the set of variogram matrices, as Bayesians do for the correlation in the 
Kriging modernization. 


Acknowledgments 

An early version of this paper was presented at the 4th Stochastic Model¬ 
ing Techniques and Data Analysis International Conference, June 1-4, 2016, 
University of Malta, Valletta, Malta, with the title Bayes and Krige: General¬ 
ities. The Authors thank both Guillaume Kon Kam King (CCA, Moncalieri) 
and Luigi Malagb (CCA, Moncalieri) for suggesting relevant references. We 
thank Emilio Musso (DISMA, Politecnico di Torino) for help with the ellip- 
topes picture. G. Pistone is supported by de Castro Statistics, Collegio Carlo 
Alberto, Montalieri, and he is a member of GNAEA-INDAM. 



12 


PISTONE AND VICARIO 


References 

[1] John Barnard, Robert McCulloch, and Xiao-Li Meng. Modeling covariance matrices in 
terms of standard deviations and correlations, with application to shrinkage. Statist. 
Sinica, 10(4):1281-1311, 2000. 

[2] Jean-Paul Chiles and Pierre Delfiner. Geostatistics. Modeling spatial uncertainty. Wiley 
Series in Probability and Statistics. John Wiley & Sons Inc., Hoboken, NJ, 2nd edition, 
2012 . 

[3] Noel A. C. Cressie. Statistics for spatial data. Wiley Series in Probability and Math¬ 
ematical Statistics: Applied Probability and Statistics. John Wiley & Sons Inc., New 
York, 1993. Revised reprint of the 1991 edition, A Wiley-Interscience Publication. 

[4] Carlo Gaetan and Xavier Guyon. Spatial statistics and modeling. Springer Series in 
Statistics. Springer, New York, 2010. Translated by Kevin Bleakley. 

[5] Tilmann Gneiting, Zoltan Sasvari, and Martin Schlather. Analogies and correspon¬ 
dences between variograms and covariance functions. Adv. in Appl. Probab., 33(3):617- 
630, 2001. 

[6] Georges Matheron. Traite de geostatistique applique. Number 14 in Mem. Bur. Rech. 
Geog. Minieres. Editions Technip, 1962. 

[7] Giovanni Pistone and Grazia Vicario. A note on semivariogram. In T. Di Battista, 
E. Moreno, and W. Racugno, editors. Topics on methodological and Applied Statistical 
Inference, Studies in Theoretical and Applied Statistics, pages 181-190. Springer, 2016. 

[8] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 
Numerical recipes: the art of scientific computing. Cambridge University Press, Cam¬ 
bridge, 1996. 

[9] Francesco Rapisarda, Damiano Brigo, and Fabio Mercurio. Parameterizing correla¬ 
tions: a geometric interpretation. IMA J. Manag. Math., 18(l):55-73, 2007. 

[10] Peter J. Rousseeuw and Geert Molenberghs. The shape of correlation matrices. Amer. 
Statist., 48(4):276-279, 1994. 

[11] Zoltan Sasvari. Positive definite and definitizable functions, volume 2 of Mathematical 
Topics. Akademie Verlag, Berlin, 1994. 

[12] Grazia Vicario, Giuseppe Craparotta, and Giovanni Pistone. Meta-models in computer 
experiments: Kriging vs artificial neural networks. Quality and Reliability Engineering 
International, 32:2055-2065, 2016. 

G. Pistone: de Castro Statistics, Collegio Carlo Alberto, Via Real Colle- 

GIO 30, 10024 Moncalieri, Italy (E-mail: giovanni.pistone@carloalberto.org). 

G. Vicario: DISMA Luigi Lagrange, Politecnico di Torino, Corso Duca degli 

Abruzzi 24, 10124 Torino, Italy (E-mail: grazia.vicario@polito.it) 



